function [ I, n ] = SimpsonIntegrationWithEps( fx, a, b, eps )
n0 = ceil((b - a) / eps^(1/4));
I1 = SimpsonIntegration(fx, a, b, n0);
I0 = SimpsonIntegration(fx, a, b, 2 * n0);
curr_eps = abs(I0 - I1) / 15;
while(curr_eps > eps)
    n0 = n0 * 2;
    I1 = I0;
    I0 = SimpsonIntegration(fx, a, b, n0 * 2);
    curr_eps = abs(I0 - I1) / 15;
end;
I = I1;
n = n0;
end

